EPJ manuscript No. 

(will be inserted by the editor) 



Monte-Carlo sampling of energy-constrained quantum 
superpositions in high-dimensional Hilbert spaces 

Frank Hantschel ^ and Boris V. Fine ^ 

Institute for Theoretical Physics, University of Heidelberg, Philosophenweg 19, 69120 Heidelberg, Germany 



o 

(N 
D 

Ph. 



26 November, 2010 

Abstract. Recent studies into the properties of quantum statistical ensembles in high-dimensional Hilbert 
spaces have encountered difficulties associated with the Monte-Carlo sampling of quantum superpositions 
constrained by the energy expectation value. A straightforward Monte-Carlo routine would enclose the 
energy constrained manifold within a larger manifold, which is easy to sample, for example, a hypercube. 
The efhciency of such a sampling routine decreases exponentially with the increase of the dimension of 
the Hilbert space, because the volume of the enclosing manifold becomes exponentially larger than the 
volume of the manifold of interest. The present paper explores the ways to optimise the above routine by 
varying the shapes of the manifolds enclosing the energy-constrained manifold. The resulting improvement 
in the sampling efficiency is about a factor of five for a 14-dimensional Hilbert space. The advantage of the 
above algorithm is that it does not compromise on the rigorous statistical nature of the sampling outcome 
and hence can be used to test other more sophisticated Monte-Carlo routines. The present attempts 
to optimise the enclosing manifolds also bring insights into the geometrical properties of the energy- 
constrained manifold itself. 
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1 Introduction 

Experimental efforts to create quantum computers come 
increasingly close to controllable manipulations of com- 
pletely isolated quantum systems consisting of the num- 
ber of g-bits of the order of 10. Although not macroscopic, 
such g-bit systems have sufficiently large Hilbert spaces, 
where it becomes increasingly difficult to generate pre- 
determined quantum superpositions. Instead, the experi- 
ments are likely to deal with the ensembles of quantum 
superpositions produced either on purpose or because of 
experimental constraints. On the theoretical side, the sta- 
tistical properties of the superpositions of quantum states 
in many-dimensional Hilbert spaces with a constraint on 
the energy value (or the expectation value of some other 
observable quantity) are also of interest for the founda- 
tions of quantum statistical physics. At issue here is the 
applicability of the Boltzmann-Gibbs statistics to com- 
pletely isolated quantum systems. 

Recently, the authors have investigated JJISJ the prop- 
erties of the so-called "quantum micro-canonical" (QMC) 
ensemble [5] of wave functions having fixed energy expecta- 
tion value (see also Refs.[4l[5l|6l[7l|8l[9]). For a Hilbert space 
of dimension TV with energy spectrum {Ei, E2, En}, 
the QMC ensemble is formally defined to include all pos- 



sible wave functions 



N 



(1) 



such that J2^=i IciP-Ei = i?av Here are the eigenstates 
of the system, Ci are the corresponding complex ampli- 
tudes, and i?av is the energy expectation value. "All pos- 
sible wave functions" in the above definition means that 
the joint probability distribution of a set of normalized val- 
ues {ci} is uniform on the manifold in the Hilbert space 
constrained by the value of i?av 

The primary objective of this paper is to explore and 
improve the numerical algorithms for sampling the QMC 
ensemble in high-dimensional Hilbert spaces. 



^ Email: F.Hantschel@thphys.uni-heidelberg.de 
Email: B.Fine@thphys.uni-heidelberg.de 



1.1 Analytical results for the QMC ensemble 

The QMC ensemble is different from the conventional micro- 
canonical ensemble: the latter limits the participating eigen- 
states to a small energy window around i?av, while the for- 
mer does not [see the discussion in Ref. 2 ]. According to 
the available analytical results[l,4,6 , the QMC ensemble 
does not lead to conventional Boltzmann-Gibbs statistics 
for either the entire isolated quantum system, or for a 
small subsystem of it. 
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For the entire quantum system with sufficiently large 
but finite Hilbert space of dimension N, the QMC ensem- 
ble typically leads to the following average occupations of 
eigenstates [UIUE] : 



N=10 



(2) 



where Pi = |cip, and A is a parameter that can be obtained 
numerically for a given energy spectrum and a given value 
of E^^. 

One can immediately see that, as Ei grows, the de- 
crease in the average occupations of the eigenstates is 
much slower than exponential. (Were it exponential, such 
a result would represent the canonical ensemble and lead 
to the conventional Boltzmann-Gibbs statistics for small 
subsystems of macroscopic systems.) 

Another remarkable feature of formula ^ is that, as 
a function of Ei, it has a pole at E\ = E^m — 1/A. As 
shown in Ref.[I], if the QMC ensemble is considered for 
a macroscopic system with the energy expectation value 
i?av equal to the average energy of the Boltzmann-Gibbs 
distribution for the same system at any experimentally 
realisable temperature, then the above pole approaches 
extremely closely to the energy of the ground state. As a 
result, the ground state acquires a macroscopically large 
occupation, which, if the QMC ensemble were realisable, 
would signify a new type of condensation phenomenon in- 
dependent of the statistics of the constitutent particles of 
the system. For quantum systems that have a large num- 
ber of states but not too large a number of constitutent 
particles, the above condensation is not sharp, but rather 
has a character of smooth crossover as a function of i?av 

For a small subsystem of a macroscopic system, the 
QMC-based result also looks unfamiliar [TJ|H]: all states of 
the subsystem except for the lowest one have the same 
occupations, while the lowest state has a higher occupa- 
tion. In other words, the subsystem appears to be in a 
weighed mixture of two conventional thermal states — the 
zero temperature state and the infinite temperature state. 

We also note that Eq. ^ already gives a good descrip- 
tion of the QMC ensemble for the values of N of the or- 
der of ten[Tl[2]. In this case, however, the deviations from 
Eq.([2]) are still noticeable for the occupations of the lowest 
(or highest) energy levels. However, these deviations can 
be accounted for with the help of the finite- corrections 
introduced in Ref.f?. 

Although the QMC-based results appear to contradict 
everyday experience, this does not mean that the QMC 
ensemble should be labelled as "unphysical" . In everyday 
experimental situations, one does not deal with completely 
isolated large quantum systems. The physical significance 
of the QMC ensemble is that it might be realisable under 
strong perturbations in isolated quantum systems with 
small numbers of particles but large numbers of quantum 
levels, for example, systems often g-bits. It is also not clear 
at present, what are the physical reasons behind the fact 
that the QMC ensemble is not realised in naturally occur- 
ring macroscopic systems. A radical possible explanation 
of this fact is based on the notion of quantum collapse. 
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Fig. 1. Unnormalised distribution of the values of iSav ob- 
tained by direct sampling of a 10-dimensional Hilbert space 
without the energy constraint. The energy levels Ei of the 
system are indicated above the x-axis. The maximum of the 
distribution is located at iJav ~ tt T^, ~ 0- 



1.2 Previous numerical investigations of the QMC 
ensemble 

The above mentioned analytical results for the QMC en- 
semble are based, in part, on approximations that need to 
be tested numerically by Monte-Carlo sampling. 

A straightforward Monte-Carlo routine involves a sam- 
pling of the entire Hilbert space with the subsequent selec- 
tion of the superpositions having the energy expectation 
values close to -Eav This routine, however, is rather inef- 
ficient. It is known from analytical calculations [HIS] that 
the occurrence of superpositions with a given energy ex- 
pectation value i?av rapidly decreases as E^^ deviates from 
the arithmetic average of all eigenenergies. As illustrated 
in Fig. [1] even in modestly large Hilbert spaces there is 
a range of values of E'av that would never appear in the 
course of such a routine implemented with realistic com- 
putational resources. The larger is the Hilbert space di- 
mension, the smaller is the fraction of all possible values 
of i?av accessible with such a routine. 

In our previous work [BIIIj in order to examine the 
statistics for any value of -Eav, we proceeded as follows. 

We note that the phases of Ci are not constrained, and 
hence have uniform probability distribution in the interval 
[0, 2tt). At the same time, the variables \ci\ can be substi- 
tuted by the eigenstate occupation variables pi = |cj;p. 
It can be shown IJ that the uniform joint probability dis- 
tribution for the normalised set of {ci} translates into a 
uniform distribution on a manifold AI defined in the space 
of variables {pi} by the following conditions: 



N 



Y,E^P^ =E, 



N 



i 

Pi > 0, Vz. 



(3) 

(4) 
(5) 



The advantage of working in the space of variables {pi} 
is that all constraints ([3][5]) are hyperplanes, and hence the 
resulting manifold has the character of hyper-polyhedron 
with flat faces. For comparison, in the space of variables Ci 
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the manifold is curved — conseqtience of the hypersphere 
normalisation condition X^f^ kiP — 1 • 

In [Tl[2] , in order to examine the statistics correspond- 
ing to any given value of i?av, we performed the Monte 
Carlo sampling within a rectangular (N — 2)-dimensional 
box around the above manifold in the subspace constrained 
by hyperplanes ([3l|4]). This sort of algorithm, while being 
able to access any value of £'av, is still relatively ineffi- 
cient. As the dimension of the Hilbert space increases, the 
volume of the enclosing box becomes exponentially larger 
than the volume of the manifold of interest, and hence the 
acceptance rate for the Monte-Carlo points becomes very 
small. The purpose of the present work was to explore how 
much one can improve the acceptance rate by reorienting 
the above rectangular box with the appropriate resizing 
of its linear dimensions, or by choosing a non-rectangular 
parallelogram-like box. 

Other algorithms to sample the QMC ensemble have 
been meanwhile proposed[6,9 and, in principle, can sam- 
ple the QMC ensemble more efficiently. These algorithms 
guide the sampling on the basis of available analytical re- 
sults. However, since some of these analytical results are 
themselves of approximate nature, it remains to be shown, 
that such algorithms lead to a fair representative sampling 
of the QMC ensemble. Given that some of the properties 
of the QMC ensemble, such as the condensation to the 
ground state, are rather non-intuitive and have unclear 
sensitivity to the numerical uncertainties, it is highly de- 
sirable to remove from the numerical investigations any 
doubts about the fair character of the Monte-Carlo sam- 
pling routines. The clear advantage of the algorithms con- 
sidered in the present paper is that, whenever they pro- 
duce sufficient statistics, this statistics is guaranteed to 
represent the true QMC ensemble. Therefore, the rela- 
tively slow algorithms described below can be used to test 
faster algorithms. 

We also note that the present effort to optimise the 
choice of the sampling box around the manifold M reveals 
interesting insights into the geometry of this manifold it- 
self. 

In a broader mathematical context, the problem of 
computing the volume of a generic convex high-dimensional 
polyhedron is known to be NP-hard|10). This and related 
optimisation problems are the subject of active ongoing 
research — see, for example, Refs. jlll[T^ll3| . In compari- 
son with the generic polyhedron problem, the case consid- 
ered in the present paper is somewhat simpler, because, 
as shown below, we can easily identify vertices, edges and 
faces of the polyhedron of interest. 

In the rest of the paper. Section \TT\ describes the basic 
idea of the Monte-Carlo sampling algorithm. Section 12.21 
describes the vertices of the manifold M, Sections 12. 31 and 
l2.4l introduce two improved algorithms. Section |3] presents 
the performance tests for the algorithms considered, and, 
finally, Section 2] summarises the results presented in this 
paper. 



2 Monte-Carlo algorithms 
2.1 Basic algorithm 

Our basic algorithm for the Monte-Carlo sampling of man- 
ifold M defined by Eqs. piS)) is based on putting a (A^-2)- 
dimensional box referred to as B around manifold M. The 
box B should lie in the the (iV— 2)-dimensional hyperplane 
A formed by the intersection of the energy and the normal- 
isation hyperplanes given by Eqs.([3]) and (0]), respectively. 
Different algorithms discussed in subsections 12.31 and 12.41 
differ only by the shape and the orientation of box B. In 
this subsection we describe all the common elements of 
these algorithms, which we call the "basic algorithm" . It 
consists of the following steps: 

1) Define the "standard" coordinate system in the space 
of variables {pi} with the origin at point (0,0, ...,0) and 
the set of TV basis vectors: ei = (1, 0, 0, 0), 62 — (0, 1, 0, .. 

ejv = (0,0,0,..., 1). 

2) Define the "modified" coordinate system, where the 
origin is chosen on one of the vertices of manifold M (that 
is, in the hyperplane A), and the basis vectors are selected 
as follows: The first two vectors are orthogonal to the hy- 
perplane A and denoted as bnrm and b^;. Vector bnmi is 
chosen perpendicular to the normalisation hyperplane 

that is, bnrm — (171,1, •■•7 1)/V^). Vector is obtained 
by orthonormalising vector (Ei, E2, E3, En) with re- 
spect to the bju-m with the help of the Gram-Schmidt pro- 
cedure. The remaining TV — 2 basis vectors, to be denoted 
as bi, b2, bAr_2, are all orthogonal to bn^n and b^;, 
that is, they lie in the hyperplane A, but they are not 
necessarily orthogonal to each other. Their specific choice 
depends on the particular realization of the algorithm as 
described in the following sections. 

3) Find all vertices of manifold AI in the standard 
coordinate system using formulas given in subsection 12.21 

4) Transform the coordinates of all vertices to the mod- 
ified coordinate system. 

5) Define the orientation of box B in the modified coor- 
dinate system; choose the linear dimensions of box B that 
are minimally sufficient to enclose all vertices of manifold 
M within the box. 

6) Randomly sample points within box B in the mod- 
ified coordinate system. 

7) Select the points belonging to manifold AI. The 
selection criterion is the following: Each sampled point 
should be transformed to the standard coordinate sys- 
tem and then accepted only if all standard coordinates 
are non-negative, meaning that the positivity condition 
^ is satisfied. 

Denoting the number of sampled points as n^, the 
number of accepted points as Ua, the (A^ — 2)-dimensional 
volume of manifold M as Vm and the volume of the sur- 
rounding box as Vb , we express the acceptance rate of the 
algorithm as follows 
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From our experience, a reasonable acceptance rate for 
practical computation should be larger than 10"''. 

Below we consider the following choices for the enclos- 
ing box B: 

1) In subsection [231 we choose the box B to be of rect- 
angular shape and optimise its orientation. Changing the 
box orientation modifies its dimensions and hence, once 
optimised, can reduce its volume. The smaller the volume 
of the box, the larger the acceptance rate r. We refer to the 
resulting algorithm as " R-algorithm" (R for rectangle). 

2) In subsection 12 .41 we enclose the manifold M within 
a non-rectangular parallelogram-like box and then opti- 
mise the choice of this box. We call the resulting algorithm 
" NR-algorithm" (NR for non- rectangular) . 

Both R- and NR-algorithms require finding the ver- 
tices of manifold M — step 3 of the basic algorithm. This 
step is described in the next subsection. 



2.2 Vertices of manifold M 

In the standard coordinate system, each vertex of the man- 
ifold M has only two non-zero coordinates fll. Moreover, 
the fact that the i-th and the j-th coordinates of the ver- 
tex are not zeroes uniquely identifies the vertex. Hence 
the vertices can be labelled by the indices of the non-zero 
coordinates. The coordinates of vertex Vij are 

= (0,...,0,K,0,...,0,Pj,0,...,0), (7) 

where 

Pi = (8) 

Pj = ^7^- (9) 

We adopt the convention that the order of indices in Vij 
is always such that Ei < Ej . Since, according to condition 
(O, both Pi and pj should be non-negative, the vertex Vij 
exists only when Ei < E^^^ < Ej. Denoting the numbers 
of energy levels below and above i?av by K and L, respec- 
tively, we can thus express the total number of vertices of 
manifold M as 

Ny^ K-L. (10) 

Obviously K + L = N. Therefore, dependent on the po- 
sition of i?av within the energy spectrum, the value of Ny 
changes between iV — 1 and (for even TV). 

In order to implement the NR-algorithm, we will also 
need to know which vertices Vk,i are connected to a given 
vertex vij by a linear edge of manifold M. The criterion 
here is that either k = i or I = j. Therefore, the total 
number Ne of edges originating from each vertex of M is 

Ne^ (K -1) + (L-l) = N ~2 (11) 

This number is thus exactly equal to the dimension of 
manifold AI for any value of i?av 




Fig. 2. (Colour online) Cartoon representing two different 
choices of boxes Bi and B2 (red rectangles) enclosing the same 
manifold M (black triangle). Obviously the volume of box Bi 
is smaller than the volume of box B2. 

2.3 Rectangular box: R-algorithm 

For the R-algorithm, we choose box i? to be a high-dimensional 
hyper- rectangle (orthotop) enclosing manifold M. In this 
algorithm, the modified basis vectors {bi, b2, bAr_2} 
specified at step 2 of the basic algorithm are chosen to be 
orthonormal, and then the edges of box B are oriented 
along these vectors. The basis vectors {bi, b2, bjv-2} 
are constructed with the help of the Gram- Schmidt pro- 
cedure, which requires N — 2 non-collinear input vectors 
gi, g2, gAr-2 defined in the standard basis. Each new in- 
put vector gi is orthogonalized first with respect to b,irm 
and b^;, and then with respect to all already available 
vectors hj. 

The orientation of the resulting dimensions of box B 
depends on the choice of the input vectors gi and on their 
sequence in the above Gram-Schmidt procedure. These 
input vectors can be chosen randomly, but such a choice 
is a priori unlikely to be optimal: indeed, it is not (see 
Section [3]). As illustrated in Fig. [2] by a two-dimensional 
cartoon, different orientations of rectangular boxes around 
a polygon can clearly lead to different box volumes, and 
hence, according to formula different acceptance rates. 
Such a difference might be more dramatic in higher dimen- 
sions. 

We have performed a partial optimisation of the box 
orientations on the basis of the following idea: As one can 
see in Fig. [51 a possibly economical way to put a box 
around the polyhedron M is to pick one face of this poly- 
hedron to coincide with one face of the box. 

To utilise this idea, we note that, according to Eqs.([31- 
[S]) each (A^— 3)-dimensional face of M is determined by the 
intersection of the hyperplane A with one of the (A^ — 1)- 
dimensional hyperplanes given by condition pi = 0, that 
is, each face can be labeled as Fi by the index of the 
corresponding eigenstate. 

Selecting the first input vector for the Gram-Schmidt 
procedure as gi = e.; causes gi to be orthogonal to Fi 
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thus guaranteeing that one of the resulting box faces will 
coincide with Fi. 

Since each face Fi is associated with a different energy 
level Ei, different faces Fi are not equivalent. Therefore, 
the question remains: Which of the N faces Fi should be 
chosen to coincide with a face of box Bl Or, equivalently: 
Which of the N natural basis vectors should be selected 
as gi? It should be further noted, that different faces Fi 
are, in general, non-orthogonal to each other. Therefore, 
once one face Fi coincides with a face of the rectangular 
box, such a coincidence cannot occur for other faces. Nev- 
ertheless, the resulting box volume may depend on the 
sequence of the remaining vectors gi used for the Gram- 
Schmidt procedure. 

Below we investigated a limited group of boxes ob- 
tained by assigning each Gram-Schmidt input vector gi 
to be equal to one of the standard basis vectors e^. The 
investigation was conducted for the following 10-level en- 
ergy spectrum representing a crudely discretised version 
of the Gaussian density of states symmetric with respect 
to zero energy: 

{Ei}i=i..io = {-0.929, -0.679, -0.466, -0.273, 

-0.09, 0.09, 0.273, 0.466, 0.679, 0.929} (12) 

It is the same spectrum as the one shown in Fig. [1] 
All possible input sequences of type 

gi = , 
g2 = , 



(13) 



lN-2 



have been tried, and the resulting box volumes for each 
sequence were obtained. All input sequences were divided 
into N groups with each group having the same standard 
basis vector Gi used as the first Gram-Schmidt input vec- 
tor gi . For each group, the minimum volume of box B was 
found over all possible combinations of vectors used in 
the rest of the Gram-Schmidt sequence. The results corre- 
sponding to i?av = — 0.3 are summarised in Table [TJ It was 
found that the global minimum volume was encountered 
in nine out of ten groups, the exception being the group 
corresponding to gi = ei with the minimum volume ex- 
ceeding the global minimum by about 20 per cent. Within 
each group, there are multiple occurrences of the mini- 
mum volume. The number of these occurrences is referred 
to in Table [T] as the "degeneracy factor" . 

For other values of E^v , the qualitative character of the 
results is the same as for i?av = —0.3. Namely, the choice 
of gi = ei for i?av < and gi = eio for i5av > does not 
allow one to reach the global minimum for the volume of 
box B. Figure [3] presents the ratio Vio/V^nin as a function 
of £^av, where Vio is the minimum volume for the group 
characterised by gi = eio and Knin is the global minimum 
volume for all groups. 

Another generic feature apparent from Table [1] is the 
greater chance of encountering the global minimum (larger 
degeneracy factor) for groups with gi = ei, such that the 



group index 

i 


min. box 
volume Vi 


degeneracy 
factor 


1 


0.0992 


6 


2 


0.0848 


1872 


3 


0.0848 


6192 


4 


0.0848 


6192 


5 


0.0848 


7056 


6 


0.0848 


7056 


7 


0.0848 


6912 


8 


0.0848 


6480 


9 


0.0848 


4320 


10 


0.0848 


72 



Table 1. Table summarising the investigation of rectangular 
box volumes for energy spectrum (|12|l with E^v = —0.3. All 
Gram-Schmidt input sequences of form (|13|) have been divided 
into groups according to the the index of the standard basis 
vector Gi used to define the first Gram-Schmidt input vector in 
gi = Bi. The left column labels these box groups. The middle 
column gives the value of the minimum box volume Vi for each 
group. The right column contains the number of times the 
minimum volume Vi is encountered within the group. 
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Fig. 3. Plot of Vio/Vm 
minimum volume for the box group characterised by gi — eio, 
and Vmin is the minimum among all boxes considered (degen- 
eracy factor). 



corresponding energy level Ei is located in the middle of 
the spectrum. 

Within each group, the influence of the choice of con- 
secutive input vectors g2,---,gAf-2 on the resulting box 
volume is complicated to describe for the present 10-level 
spectrum, and we do not attempt it here. We have, how- 
ever, also conducted similar investigations for smaller spec- 
tra with N = 5 and = 6, where a clear picture emerged. 
In the both N = 5 and A^ = 6 cases, we have found that 
(i) the last input vector gAf-2 does not confine the volume 
in any way; and (ii) the minimal box volume is realised, 
whenever neither of the vectors gi, gNs is assigned to 
be ei or e^. Therefore, for each input group where the 
global minimum volume appears (that is, for gi — Gi^ 
with ii — 2, ...,N — 1), the degeneracy factor is equal to 
3(A^ — 3)!. While such a result is valid for A^ = 5 and 
A^ = 6, it does not hold for A^ = 10 (see Tabled]). 

Our final prescription for the R-algorithm is the fol- 
lowing: At step 2 of the basic algorithm, we construct 
orthonormal modified basis vectors bi, bjv-2 using the 
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Fig. 4. (Colour online) Cartoon of parallelogram-like box B 
(red parallelogram) enclosing the manifold M (black triangle). 
The large dot denotes the "origin vertex" of M, while the ar- 
rows indicate the basis vectors used to create B. 

Gram-Schmidt orthogonalization procedure with input vec- 
tors gi = e^j^,..., gN-2 = ei„_2, where none of the vectors 
e^j. are equal to either ei or e^v, and otherwise vectors e^j. 
appear in the order of the proximity of the correspond- 
ing energies Ei^ to the arithmetic average of all energies 
in the spectrum. At step 5, we choose a rectangular box 
with the edges oriented along vectors bi, ...,bjv_2- 

In each calculation, the above prescription can be par- 
tially controlled by doing random sampling of a large num- 
ber of Gram-Schmidt input sequences {gfc = 6^^.} and 
then checking that the volumes of the resulting boxes are 
not smaller than the volume of the "prescription box" . 
All such tests performed in the specific cases presented in 
Section [3] have supported the above prescription. 




Fig. 5. Volumes of parallelogram-like boxes Vb for different 
choices of the origin vertex Vij for the Gaussian spectrum of 
A'^ = 10 and i?av = 0. Coordinates i and j are the indices of 
the chosen vertex Vij. choice. 



-Eav 


K 


L 


Vmin 


Rvol 





5 


5 


9.93 10" 




3.11 


-0.04 


5 


5 


9.86 10" 


3 


3.11 


-0.08 


5 


5 


9.6410" 


3 


3.11 


-0.10 


4 


6 


9.2710" 


3 


2.97 


-0.46 


3 


7 


8.20 10" 


3 


2.55 


-0.5 


2 


8 


4.60 10" 


3 


1.88 



Table 2. Properties of the volumes of parallelogram-like boxes 
B for the 10-level energy spectrum H12|) as a function of -Bav- 
The parameters of the table are defined in the text. 



2.4 Non-rectangular box: NR-algorithm 

The NR-algorithm uses non-rectangular parallelogram-like 
box for B. In order to define such a box, one needs to spec- 
ify one of the box vertices, which we refer to as the "origin 
vertex" , and the N — 2 non-orthogonal edges originating 
from this vertex. As indicated by Eq.([TO]), each vertex of 
manifold M also has exactly N — 2 edges originating from 
it. Therefore, the origin vertex of box B is chosen to coin- 
cide with one of the vertices of manifold M, and then the 
corresponding N — 2 edges of the manifold M determine 
the edge directions for box B. A two-dimensional cartoon 
of such an arrangement is shown in Fig. [J] The origin 
vertex also becomes the origin of the modified coordinate 
system, and the vectors of the modified basis are oriented 
along the same TV — 2 edge directions of box B. 

Such a general algorithm still leaves the freedom of 
choosing a vertex of manifold M as the origin vertex of 
box B. The number of vertices of manifold M is given 
by Eq.(I5]). Below we explore the dependence of the box 
volume on the choice of the origin vertex for the energy 
spectrum A possible algorithm for the calculation 

of the volume of a parallelogram-like box is given in Ap- 
pendix 

As explained in Section [^21 each vertex is labeled as 
Vij, where the two indices are such that Ei < Es,-v and 
Ej > Eav ■ Figure [5] presents the results of box volume 
calculations for for _Eav = 0. In this case, E^ < i?av < E^. 



Therefore, there exist 25 vertices Vij with index i tak- 
ing one of the values {1, 2, 3, 4, 5}, while index j may take 
values {6, 7, 8, 9, 10}. The results are presented as a three- 
dimensional plot, where two horizontal axes represent ver- 
tex indices i and j and the vertical coordinate represents 
the resulting volume. In this figure, the ratio i?voi of the 
maximum volume Vmax over the minimum volume Knin is 

= = 3.10627 (14) 

As one can see from Fig. [5l there exist nine cases realizing 
the minimum volume. They correspond to vertices of the 
form vi,j with arbitrary j, or Wi,io with arbitrary i. The 
maximum volume corresponds to only one vertex v^^e with 
two respective levels bracketing i?av from below and above. 
Qualitative arguments explaining, why it is expected that 
such a situation is special, are given in Appendix |B] 

We have further investigated the characteristics of box 
volumes, for several values of i?av and collected the results 
for Vmin and i?voi in Table [51 This table demonstrates that 
(i) the ratio i?voi is the same for different values of E^v but 
agreeing values of K and L; and (ii) i?voi decreases with 
decreasing number of vertices (equal to KL). 

Our final prescription for the NR-algorithm is to choose 
vertex vi^n as the origin vertex. (Since, from the compu- 
tational viewpoint, the number of possible choices of the 
origin vertex is not large, one can always check the above 
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prescription by computing the box volumes for all possible 
choices.) 



3 Performance of R- and NR-algorithms 

This section compares the performance of the algorithms 
presented in section [^751 and by looking at their accep- 
tance rates. The following conditions were chosen: 

1) Number of states N was equal to 10, 12, or 14 ; 

2) Energy spectra represented discretised realisations 
of the Gaussian density of states similar to spectrum (|12l) 
with arithmetic average of all energies equal to zero and 
root- mean-squared deviation from zero equal to 1/ . 

3) The average energy £^av was chosen between and 
O.SEi. {El is the minimum eigenenergy of the spectrum.) 

Table [3] contains the acceptance rates for three differ- 
ent algorithms: (i) non-optimised R-algorithm with ran- 
dom isotropic selection of Gram- Schmidt input vectors 
gi, ...,gAr_2; (ii) optimised R-algorithm with prescription 
for the input Gram-Schmidt vectors given at the end of 
subsection l2.3l and (iii) optimised NR-algorithm with pre- 
scription for choosing the origin vertex given at the end 
of subsection [231 We use non-optimised R-algorithm as a 
benchmark for inefficient choice of box B, which allows us 
to judge the effect of optimising box B. 

As one can see from Table [31 the value of acceptance 
rates r for all algorithms decreases rapidly with the in- 
crease of to the extent that the sampling of Hilbert 
spaces with dimensions larger than 14 appears impracti- 
cal with the present algorithms. 

For a given spectrum, the acceptance rates also de- 
crease with the increase of the ratio E^x/Ei. This fact is 
further illustrated in Fig. [HI for the R-algorithm applied to 
the spectrum with N — 10. An interesting detail revealed 
by Fig. [6l is that the acceptance rate does not depend on 
E^^ for El < E^av < E2. 

Comparing the relative performance of the three algo- 
rithms, it can be observed from Table [31 that the optimi- 
sation of R- or NR- algorithm does not result in significant 
increase of the acceptance rates for the values of i?av close 
to the centre of the spectrum. However, as the average 
energy deviates from the centre the optimisation leads to 
the increase of the acceptance rates by about a factor of 5. 
The acceptance rates for optimised R- and NR-algorithms 
are close to each other for all combinations of parameters 
considered. 

The decrease of the acceptance rates with the depar- 
ture of i?av from the centre of the spectrum, as well as the 
inscnsitivity of the algorithms to the optimisation for i^av 
close to zero, can probably be attributed to the fact that, 
according to section 12.21 the number of vertices of mani- 
fold M is maximal — equal to A^^/4 — for iSav = 0, and 
then it quickly decreases towards A^ — 1 as ii^av deviates 
from zero. Larger number of vertices, supposedly indicates 
that manifold M has more even shape, which fills a larger 
volume fraction of any box around it irrespective of the 
box orientation. On the contrary, the smaller number of 
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Fig. 6. Plot of the acceptance rate r of the R-algorithm versus 
the average energy iSav for the 10- level energy spectrum p2p . 

vertices may imply more uneven shape of manifold M oc- 
cupying a smaller fraction of any reasonably-shaped box, 
with the occupied fraction being strongly dependent on 
the box orientation. 



4 Conclusions 

This paper described and optimised two algorithms for 
performing Monte Carlo sampling of quantum superposi- 
tions in high-dimensional Hilbert spaces under the con- 
straint on the energy expectation value. The two algo- 
rithms are distinguished by the shape of Monte- Carlo sam- 
pling boxes: either rectangular or parallelogram-like. The 
optimisation included finding the box orientations allow- 
ing smaller box sizes and, therefore, larger Monte-Carlo 
acceptance rates. Both algorithms were found to exhibit 
similar performances in the optimised form. The bene- 
fit of the optimisation in terms of algorithm's acceptance 
rate was found to be small for the values of E^v close to 
the centre of the Gaussian-like spectra considered but in- 
creased by about factor of 5 as E^v moved to the wings of 
the spectra. We have found that the largest dimension of 
Hilbert space that can realistically be explored with these 
algorithms is 14. 

Our efforts to optimise the choice of the sampling boxes 
also indicate that for E^v close to the centres of energy 
spectra, the manifolds of interest are more even-shaped, 
while, in the more interesting regime of i^av being on the 
wings of the energy spectra, the resulting manifolds are 
more uneven and difficult to sample. 



A Volume of high-dimensional 
parallelogram-like box 

In order to determine the volume of an (A^— 2)-dimensional 
parallelogram-like box, the following algorithm can be ap- 
plied iteratively. Each iterative step of this procedure is 
illustrated by Fig. [3 

The numerical routine begins by selecting one vertex 
as the origin and then determining all edges connecting 
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Table 3. Acceptance rates r of R- and NR- algorithms for several spectra and several values of iJav as described in the text. 
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Fig. 7. Cartoon of a two-dimensional parallelogram illus- 
trating the calculation of the volume of high-dimensional 
parallelogram-like box. The volume (area) of this parallelo- 
gram is 62 /ii, where hi is the height and 62 is the base side as 
indicated. 

the origin vertex to the N — 2 neighbouring vertices. The 
vector basis which is defined by the set of all edges is 
denoted as {bi, b2 . . . bjv_2}. 

The first iteration takes the subset of the basis vectors 
{b2 . . .b7v-2} and finds a normalised vector Si which is 
orthogonal to each vector of the above subset. This is done 
by solving a system of linear equations for the orthogonali- 
sation conditions. Once Si is obtained, the following height 
parameter can be defined: 

/ii = (bi-si) (15) 

The second iteration finds a normalised vector S2 which 
is orthogonal to {bs . . . bAr_2! Si} and then finds the cor- 
responding /l2. 

The iterations are continued until /iat-s is found. The 
resulting volume of the box is then given by 

VB-|bAr_2|- (16) 

B Maximum volume box in the NR-algorithm 

According to the results of Section l^?^ the maximum sam- 
pling box volume for the NR-algorithm corresponds to the 
choice of the origin vertex Vi_j such that Ei and Ej bracket 
the average energy E^v- This choice is special because of 
the geometrical factors described below . 



The vertices connected by the common edges to Vij 
form two groups: the first group consists of vertices Vi^m 
(the same first index as for the origin vertex), while the 
second group consists of Vnj (the same second index as for 
the origin vertex). As discussed in Section [221 each vertex 
Vi^m in the first group has only two non-zero coordinates 
Pi and Pm, where: 

= (17) 

Pm = fcf ■ (18) 

In this group, the typical situation is that Ei is much closer 
to £'av than Em- Therefore, pi is close to 1, and Pm is close 
to zero. In other words, most vertices of the first group 
cluster around the point (0, 0, 1^, 0, 0), where the no- 
tation li indicates that pi = 1. Likewise the vertices of 
the second group cluster around (0, 0, 1^, 0, 0). This 
clustering results in relatively small angles between all ba- 
sis vectors connecting the origin vertex to the vertices of 
the same group, which presumably leads to a non-optimal 
choice of box B around AI. 
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